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Abstract. Firstly, we give a historical overview of attempts to incorporate magnetic fields into 
the Smoothed Particle Ifydrodynamics method by solving the equations of Magnetohydrody- 
namics (MHD), leading an honest assessment of the current state-of-the-art in terms of the 
limitations to performing realistic calculations of the star formation process. Secondly, we dis- 
cuss the results of a recent comparison we have performed on simulations of driven, supersonic 
turbulence with SPH and Eulerian techniques. Finally we present some new results on the re- 
lationship between the density variance and the Mach number in supersonic turbulent flows, 
finding crf^p ~ ln(l + b^A^^) w ith h = 0.33 up to Mach 20, consistent with other numerical 
results at lower Mach number l|Lemaster fc Sto"n3 l2008h but inconsistent with observational 
constraints on cTp and JVi in Taurus and IC5146. 
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1. Introduction 

Magnetic fields and turbulence are thought to be two of the most important ingredients 
in the star formation process, so it is critical that we are able to model the effect of 
both in sta r formation simulations. Sm oothed Particle Hydrodynamics (SPH, for recent 
reviews see lMonaghanll2005t |Pricell2004[) . since it is a Lagrangian method where resolution 



follows mass, is a very natural method to use in order to model star formation. However 
the introduction of magnetic fields in SPH has a somewhat troubled history, and there 
are perceptions that the explicit use of artificial viscosity terms in order to treat shocks 
is too crude to enable accurate modelling of supersonic turbulence. We will discuss both 
of these aspects in this talk. 



2. Magnetic fields in SPH 

2.1. A prehistory of MHD in SPH 
Magnetic fields were introduced in one of the founding SPH papers bv lGingold fc Monaghan 



(1977.) , though a detailed investigation did not fo llow until around a decade later with 
the publication of Graham Phillips' PhD work (jPhillips fc Monaghanl Il985l ) (PM85), 



with successful application of the method to star formation problems involving the col- 
lapse of isothermal, non-rota ting, magnetised clouds (e.g. Phillips 1986j3. Indeed, it was 
Phillips fc Monaghanl (|l985[) who first discovered that the equations of motion for MHD 



f Graham Phillips now pursues a rather different career as presenter of the science show 
'Catalyst' on Australian television (I suspect this was easier than getting MHD in SPH to 
work) . 
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- when expressed in a momentum-conserving form - were unstable to a clumping or 'ten- 
sile' instability whereby particles attract each other unstoppably along field lines, that 
occurs in the regime where the magnetic pressure exceeds the gas pressure (i.e., /3 < 1). 

Whilst Phillips & Monaghan proposed a way of achieving stability, another decade 
followed before the matter was investigated in any i nore detail, b y another PhD student 
of Joe Monaghan's, Joe Morris. Morris' PhD thesis (" MorrislflQQGl) contained detailed sta- 
bility analysis of the equations of the now-named Smoothed Particle Magnetohydrody- 
namics (SPMHD) and proposed a way of avoiding the instability by using more accurate 
derivatives for the anisotropic part of the magnetic force in place of exact conservation of 
momentum. Around the same time iMe glickil (|1995I ) applied a non-conservative SPMHD 
formulation to magnetic fields in the Galaxy. The period after these works saw a few but 
limited applications of SPMHD to 'real' astrophy sical problems, including notably a n 



application to magnetic fields in galaxy cluste rs by Dolag. Bartelma nn & Lesch ( 1999^ 3 



and at least one foray into star formation bv iByleveld k. Pongraci c (1996). However it 
would be fair to say that there was no real consensus on a standard approach and that 
many problems remained. 

2.2. 21st century Smoothed Particle Magnetohydrodynamics 
Borve. Omang Sz TrulsenI (l200lh suggested an alternative to Morris' stabilisation based 



on explicitly subtracting the non-zero source term (-BV • B) that is the cause of the 
stability problems in the momentum-conserving formulation. They also proposed a regu- 
larisation scheme for the SPH particle distribution which, though complicated, was found 
to give ex cellent results on MHD sho ck problems by the focussing of resolution where it 
is needed. IPrice fc MonaghanI ( 2004a ) adopted (but later abandoned) a stability fix pro- 



posed bv iMonagha n (200Q:) and showed how the dissipative terms could be formulated 
for treating MHD shocks and other di scontinuities in a standard SPMHD approach (that 
is, without regularisation) . Paper H (jPrice fc Monag halll2004hh showed how the terms 



necessary for strict conservation in the presence of a spatially variable smoothing length 
could be included by deriving the SPMH D equations of motion from a variational princi- 
ple. At the same time lBorve et al. (2004) undertook a detailed stability analysis showing 



that indeed their source-term approach was indeed stable in more than one dimension. 

With these improvements meant that SPMHD could be succ essfully run on many of th e 
test problems used to test grid-based MHD codes (Paper HI. IPrice fc Monaghan 2005 ). 



However, accuracy in realistic problems was found to depend on the accuracy with which 
the divergence-free condition on the magnetic field - perhaps the key problem for any 
numerical MHD scheme - since the V ■ B = constraint enters the ideal MHD equations 
only as an initial condition: 

— = Vx(vxB); _(V.B)^0. (2.1) 

As aptly summarised by Mike Norman (this proceedings) , approaches to the divergence 
constraint can be roughly categorised into "ignore", "clean" or "prevent". Many test 
problems for MHD can be run quite happily (in SPH or otherwise) with the "ignore" 
approach: simply evolve the induction equation (|2.ip and cross your fingers (I refer to this 
as the "hope and pray" method). For realistic applications this is almost never a good idea 
and in SPMHD simulations of s tar formation manifests as 'explodin g stars' due to large 
errors in the magnetic force fsee lPrice fc Federrath l2009l for a figure) . iPrice fc MonaghanI 



(,2005. ) examined a range of cleaning methods for the divergence constra int, including 
most promisingly the hyperbolic/parabolic cleaning scheme proposed bv iDedner et al] 

f Harald Lesch is now also a science show presenter, on German television. 
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( 20021) . However this was not found to be sufficient to prevent divergence errors generated 
in realistic flows, least of all during protostellar collapse. 

2.3. The Euler Potentials 

Thus, Price fc Batd (j2007l) and lRosswog fc Pried pOOTI ) adopted a preventat ive approach , 
based on writing the magnetic field in terms of the 'Euler Potentials' ([Stern Il976l) . 
variously known as 'Clebsch' or 'Sakurai' variables: 



B = Va X V/3, 



(2.2) 



with which the divergence constraint is satisfied by construction. Furthermore the ideal 
MHD induction equation takes the particularly simple form 



dt 



(2.3) 



corresponding physically to the advection of magnetic field lines by Lagrangian particles. 
This is then combined with a standard SPMHD approach by computing B according to 
Eg. 12.21 and using the standard equation of motion with either the Morris or iBorve et al. 
(|200ll) force stabilisation (thus is is usually trivial to revert to a B-based formulation - 
that is, solving Eq. 12.11 instead of Eqs. I2.2H2.3I - for comparison purposes. 

Adoption of this approach meant that we were able to study realistic star formation 
problems in 3D, including the m agnetised collapse of a rotating core to form single and 
binary stars (IPrice &: Batdl2007l ) and the effect of rnagne tic fields on star cluster forma- 
tion from turbulent molecular clouds (IPrice fc Batd[2008h . most recently including both 
radiation and magnetic fields ( Price fc Batdl2009f) . Around this time magnetic fields w ere 
also implemented in the widely used SPH code GADGET (iDolae fc Stasvszvnll2009l) . 

From these we have learnt that magnetic fields can significantly change star forma- 
tion even with weak (supercritical ) fields, s ignificantly affecting the formation of cir- 
cumstellar discs a nd binary stars ( Price fc Bate 2007„ similar to the find i ngs of many 



other authors, e.g. iHennebelle fc Tevssieill2008t iHennebelle fc Ciardill2009t iMellon fc fj 



2008t iMa chida ei^ al.' '2008') and having a profound effect on the star formation rate 



( Price fc Bate. 2008, ,2009.1 

However there are important limitations to the Euler Potentials (EPs) formulation. 
The most important are that: 

(a) They can only be used to represent topologically trivial fields, 

(6) As a consequence of a) , the wind-up of a magnetic field in a rotating or turbulent 
flow cannot be followed indefinitely, 

(c) It is difficult to incorporate single fluid non-ideal MHD effects such as resistivity. 
The topological restrictions are a consequence of the fact that the helicity A • B is identi- 
cally zero in the EP representation. The lack of winding processes is easily understood as 
a consequence of Eq. 12.31 Since the potentials are simply advected by the particles, the 
reconstruction of the magnetic field via Eq. 12.21 relies on a one-to-one mapping between 
the initial and final particle positions. Thus, for example, a complete rotation of an in- 
ner ring of particles will return the magnetic field to its original configuration, whereas 
physically the field should 'remember' the winding [in the talk this was illustrated by a 
dance with J. Monaghan which is difficult to insert in the proceedings]. 

In practice this means that in using the EPs one is limited to the study of problems 
where the magnetic field is initially simple (e.g. imposed as a constant field) and at some 
scale the field winding will be lost (underestimated) in the calculations. The difficulty in 
representing physical resistive processes is also an issue given the importance of non-ideal 
MHD effects for star formation. 
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2.4. Why don't you just use the Vector Potential? 

At one of the last major computational star formation meetings, held at the Kavli In- 
stitute for Theoretical Physics during 2007, Axel Brandenburg suggested that the use of 
the vector potential, with the gauge set to give Galilean invariance, may provide similar 
advantages in terms of the divergence constraint without th e topologic al restrictions of 
the Euler potentials. Three years and 25 pages of pain later (|Pricell2010l ) we had derived 



the Uhimate Vector Potential Formulation'^^ for SPMHD. It was beautiful, derived el- 
egantly from a variational principle, the method was exactly conservative, with a brand 
new force equation that in principle might not suffer from instabilities since the diver- 
gence constraint was built-in... could it be that we had uncovered the ultimate method for 
MHD in SPH? Unfortunately the vector potential approach failed spectacularly in prac- 
tice due to two problems: 1) The same old force-related instability in 2D and 3D, which 
could, however, be cured by reverting to the stable force formulations discussed above, 
though with more difficulty than in the standard SPMHD case, and 2) An instability 
peculiar to the 3D vector potential related to the unconstrained growth of non-physical 
components of A, occurring, for example, when a 2D test problem like the Orszag-Tang 
vortex is run in 3D. Further investigation suggests the latter is related to not enforcing 
the gauge condition V • A = on the vector potential - that is, simply pushing the 
divergence problem one level up. Thus, we must unfortunately report that the answer to 
the vector potential question is that it not a good approach for SPMHD. 

2.5. Alternatives to the Euler Potentials 

The most promising approach to lifting the restrictions of the Euler Potentials whilst 
preserving the benefits involve moving to a generalised formulation, where the most 
general involves three sets of potentials, in the form 

B = Vai X V/3i -I- Vq!2 x V/32 + Vas x V^3, (2.4) 

essentially taking advantage of the fact that, whilst sets of potentials cannot be added 
together, the fields can. Thus one can represent arbitrary field geometries. The reason 
for the three sets is because one of the potentials (3 is essentially preserving one of the 
initial spatial coordinates so that a Jacobian map can be constructed - for example 2D 
fields are represented by a = a{x, y) and /3 — zq. Thus in the most general case all three 
coordinates of the particles (/3i, /32, /Sa) = {xo,yo,ZQ) are used to compute the mapping. 
Better still, adopting this more general formulation means that the potentials can be 
re-mapped to a new set whenever the gradients break down, meaning that physical field 
winding can be followed indefinitely, and with explicit control on the amount of numerical 
dissipation. Watch this space. 

An alternative is to return to the 'clean' approach , desp ite the only limited success 
of the formulations considered bv lPrice fc Monaghanl (2005), at least for star formation 



problems. Use of cleaning methods has found more success in cosmological applications 



leaning 

(|Dolag fc Stasvszynl 20091 ) and there are further possibilities fo r improving the accuracy 



of the projection methods discussed in I Price fc Monaghanl ([20050 that offer some promise 



3. Supersonic turbulence with SPH 

Progress on solving the issues with magnetic fields in SPH generally slows due to the 
need to periodically address concerns over the core SPH method. Most recently this has 
involved a rather infiamed and unhelpful discussion of Kelvin-Helmholtz instabilities in 
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Figure 1. Time- averaged kinetic-energy power spectra (left), and the time evolution of the 
maximum dens ity (right) at three differe nt resolutions in the SPff vs. grid driven turbulence 
comparison by iPrice fc Federrattil (|2010f ). Comparable results on volume- weighted quantities 
such as the kinetic energy power spectrum requires roughly Uceiia ~ nparUcUs (left), but SPH 
resolves maximum densities at 128^ similar to a grid at 512^ (right). 



and which many authors are now using as the primary justification for the often 
fruitless pursuit of all manner of alternative particle methods. Computational star for- 
mation is no stranger to such debates, though thankfull y at least the heated discussions 
surrounding 'artificial fragmentation' ( Klein et al.ll2004 ) went quietly away, leading to a 



more helpful consensus on the implications of the physical approximations (e.g. use of 
a barotropic equation of state vs. radiation hydrodynamics) rather than arguments over 
numerics. However there remain sharp disagreements over num erics in some areas, for ex- 



ample the follo wing in a discussion o f turbulence simulations by lBallesteros-Paredes et al 
(|2006l) given in lPadoan etail (|2007t ): 



"SPH simulations of large scale star formation fields to date fail in all three 
fronts: numerical diffusivity, numerical resolution and presence of magnetic fields. 
This should case serious doubts on the value of comparing predicti ons based on 
SPH simulations with observational da ta (see also 'Agertz e t aLll2007() ". 
(note the completely irrelevant reference to iAgertz et al. 200 7). Such concerns over turbu- 
lence helped motivate at least two major comparison projects - the 'Potsdam' comparison 
(|Kitsionas et al.ll2009l ) and the 2007 KITP comparison project, both studying the decay 



of supersonic turbulence from pre-evolved initial conditions. 

Given that comparisons of decaying turbulence are necessarily limited to a few snap- 
shots and that there are issues surrounding the setup of initial conditions between codes, 
we additionally undertook a side-by-side comparison of driven turbulence, starting from 
simple initial conditions and over many crossing tim es, using just two codes - D. Price's 
SPH code PHANTOM and the grid-based flash code ( Price fc Federrathll2010l ). The dev- 



astatingly obvious conclusion from such efforts was that to achieve comparable results 
between codes required ~ wait for it - comparable resolutions. Actually the point was 
rather more subtle, because "comparable resolution" depended on what quantity one 
was interested in. Thus, for example to obtain comparable power spectra in a volume- 
weighted quantity such as the kinetic energy required roughly equal numbers of SPH 
particles to Eulerian grid cells (and was thus about an order of magnitude more costly 

f Following the publication of lAgertz et al.l (120071 ). The issue itself has nothing to do with 
Kelvin-Helmholtz instabilities but rather to the (non-)treatment of contact discontinuities in 
standard SPH schemes, in this case manifested when a KH instability problem is run ac ross a 
contact discontinuity. Adding a small amount of artificial conductivity is a simple fix (see lPric"3 
I2008f ). 
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Figure 2. The measured relationship (points) between the standard deviation in the log density, 
Ga and the RMS Mach number, JVi, from a series of calculations of solenoidally-driven, supersonic 
turbulence at various Mach numbers and at different resolutions in both SPH particles and on 
the grid to which the density field is interpolated. Overplotted is the standard relationship, 
Eg. 14.21 with b = 0.33 (lower) and b = 0.5 (above). Also shown is the best fit obtained by 
iLemaster fc Ston3 ([2008 ) at lower Mach numtser, Eq. 14.31 which differs only slightly from the 
standard relation with 6 = 0.33. Error bars show the la deviations from the time-averaged 
values. 



for SPH in terms of cpu time). On the other hand, the maximum density achieved in 
the SPH calculations at 128'^ particles was roughly equivalent to that reached in the 
grid-based calculations at 512^ (thus being about an order of magnitude more costly 
for flash). Quantities involving a mix of the density and velocity fields - such as the 
spectra of the mass weighted kinetic energy p^^^v - showed intermediate results. Other 
quantities such as structure function slopes were n ot well converged in eithe r code at 
512'^ computational elements. We refer the reader to [Price fc FederrattJ (|201(J) for more 
details. 



4. The density variance Mach number relation in supersonic 
turbulence 

One of the questions to arise from the turbulence comparison project was the ques- 
tion of exactly how the width of the density Probability Distribution Function (PDF) 
in driven, supersoni c turbulence depends on the Mach number. Early calculations by 
Padoan et al. ( 1997[ ) found/assumed a relationship between the linear density variance 



(Tp and the Mach number of the form 

al = b^M^, (4.1) 
such that in the log-normal distribution the standard deviation is given by 

<j'j ^ln{l + b^M^), (4.2) 
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Figure 3. The directly measured (lower points) and inferred (assuming a log-normal distribu- 
tion: a'p — exp((T^) — 1), higher points) relationship between the standard deviation in linear 
density, o-p, and the RMS Mach num ber, M, in su personic turbulence simulations, together 
with the best-fitting relationship from iLemaster fc S tone (2008) ( solid li n e, sim ilar to & = 0.33 
in E g. 14.111 and th e observational results from IC5146 bv iPadoan et al.l (j 19971 ) and in Taurus 
from lBruntI (|2010l '). 



where Cs = crx^p is the standard deviation in t he log of the density and h is a par ameter 
that was determined empiricaUy to be 6 '--^ 0.5 (jPadoan. Nordlund fc Joneslll997l ). Later 
authors have found considerable variation in the measured value of b based on indi vidual 
calculations performed at a fixed Mach number. Recently iFederrath et al. I (|2010l) have 
also shown that the value of 6 also depends strongly on particulars of the fourier-space 
driving, specifically the balance between energy in solenoidal and compressible modes 
(with b K, 0.33 for purely solenoidal and fowl for purely compressible driving), which 
they relate physically to whether or not turbulence is driven by compressions or shear 
flows. 

However, relatively few studies have attempted to constrain the relation over a range 
in Mach number, though the early e mpirical (but unpu b hshed ) work by Padoan et al. 
was calibrated by such calculations. iLemaster fc Ston3 (j2008l ) have recently done so, 
though only in a limited range of Mach numbers (up to 6) , deriving a best-fitting 

relationship 

erf = -0.72 In (l-h0.5>l2) 0.20. (4.3) 

However we have re cently developed a technique for inferring the true 3D variance dp from 
observational data ( Brunt et al.l 2O10l) that, together with the m easured Ma ch number, 
indicate 6 « 0.5±0.05 in the Taurus Molecular Cloud at M ach 20 (|B runt"2010'), similar to 
an earlier result derived at in IC5146 at Mach 10 by P adoan. Jones k No rdlund (19971). 
In order to make comparisons with observations it is therefore necessary to pin down the 
relationship up to these Mach numbers. 

We have thus performed a series of calculations with RMS Mach numbers in the range 
1 — 20, the results of which are shown in Fig. [2] in terms of CTs, measured directly as the 
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(volume-weighted) standard deviation in the (gridded) log density field with no assump- 
tions regarding log-normality or otherwise. The difficulty in comparing with observations 
is that the variance in the linear density field is severely truncated by finite resolution 
effects, but it is this quantity that is measurable observationally (though it will also be 
underestimated). One way of estimating the true value is to simply assume that the PDF 
is log-normal, i.e., = exp((T^) — 1. The calculations were observed to approach this 
relationship as the resolution - either in terms of the number of particles or the grid to 
which they are interpolated - is increased. The standard deviation in the linear density 
is what is measured in the observational data, a comparison with which is shown in 
Fig. El Whilst we find a good fit to Eq. 1121 with b = 0.33 (and also to the LS08 best fit, 
Eq. 14. 3p for the purely solenoidal driving we have employed, it is clear that the density 
variances measured observationally lie significantly above those in the calculations, with 
b « 0.5 (and this noting that our a ssumption of log - norm ality will overestimate ap, if 
anything). In the heuristic model bv lFederrath et ah ( 20081 ) this can be explained if the 
is a significant compressive component to the driving mechanism, or alternatively due to 
the fact that star-forming molecular clouds are quite obviously self-gravitating which is 
not accounted for in pure turbulence models. 



5. Conclusions 

In summary: 

• Being a television presenter is easier than getting MHD in SPH to work 

• Development of MHD in SPH would proceed faster in the absence of arguments over 
basic misunderstandings of SPH which require lengthy comparison projects to resolve. 

• Simulations with supercritical field strengths already tell us that magnetic fields can 
significantly change star formation even with weak fields. 

• SPH and grid codes agree very well on the statistics of supersonic turbulence pro- 
vided comparable resolutions are used (for power spectra this means riceiis ~ ^particles , 
but SPH was found to be much better at resolving dense structures). 

• We have constrained the density variance - Mach number relation in supersonic 
turbulence to = ln(l -I- b^M^), with b « 0.33 up to Mach 20, but observed density 
variances suggest b ~ 0.5, higher than can be produced with purely solenoidally-driven 
turbulence - rather, some form of compressive driving or additional physics such as 
gravity is required to explain the observations. 
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